
/* 
* sepprfrc.c - This file is a part of the sep-library 
*
* Copyright (C) 2008 Jesper Schmidt Hansen 
* 
* License: GPL - see COPYING for copying conditions.
* There is ABSOLUTELY NO WARRANTY, not even for MERCHANTIBILITY or
* FITNESS FOR A PARTICULAR PURPOSE.
*
* Contact: schmidt@zigzak.net
*/


#include "sepprfrc.h"
#include "sepinit.h"



void sep_force_pair_brute(seppart *ptr, const char *types, double cf,
			  double (*fun)(double, char), sepsys *sys, 
			  sepret *retval, const int opt){
  int n,m,k;
  double r[3], r2, ft, f[3], cf2; 
  
  cf2 = cf*cf;

  for (n=0; n<sys->npart-1; n++){
		
    int i = ptr[n].molindex;
					
    for (m=n+1; m<sys->npart; m++){

      if ( opt == SEP_NEIGHB_EXCL_BONDED && sep_bonded_direct(ptr, n, m) == 1 ){
	continue;
      }
      else if ( opt == SEP_NEIGHB_EXCL_SAME_MOL && 
		ptr[n].molindex == ptr[m].molindex &&
		ptr[n].molindex != -1 )
	continue;


      if ( (ptr[n].type == types[0] && ptr[m].type == types[1]) || 
	   (ptr[n].type == types[1] && ptr[m].type == types[0]) ){
				
	r2 = 0.0;
	for (k=0; k<3; k++){
	  r[k]  = ptr[n].x[k]-ptr[m].x[k];
	  sep_Wrap( r[k], sys->length[k] );
	  r2   += r[k]*r[k];
	}
	if (r2 < cf2){ 

	  // Force between particles	
	  ft = (*fun)(r2, 'f'); 
					
	  for (k=0; k<3; k++){
	    f[k] = ft*r[k];
						
	    ptr[n].f[k] += f[k];
	    ptr[m].f[k] -= f[k];
	  }
				
	  // Energy 
	  retval->epot += (*fun)(r2, 'u');
					
	  // Config. part of the stress/pressure tensor 
	  for (k=0; k<3; k++)
	    for ( int kk=0; kk<3; kk++ ) retval->pot_P[k][kk] += f[k]*r[kk];
	
	  // Force between molecules
	  if ( sys->molptr->flag_Fij == 1 ){
	    int j = ptr[m].molindex;
	    if ( i != -1 && j != -1 ){
						
	      for ( k=0; k<3;k++ ){
		// pressure tensor	
		sys->molptr->Fij[i][j][k] += f[k];
		sys->molptr->Fij[j][i][k] -= f[k];
		// Couple tensor
		sys->molptr->Fiajb[n][j][k] += f[k];
		sys->molptr->Fiajb[m][i][k] -= f[k];	
	      }
	    }	
	  }
					
	} // End of if ( r2 < rcut )
      }
    }
  }

}
	

void sep_force_pair_neighb(seppart *ptr, const char *types, double cf,
			   double (*fun)(double, char), sepsys *sys, 
			   sepret *retval) {
  int i1, i2, n, k;
  double r2, ft, f[3], r[3];
  const double cf2 = cf*cf;

  for (i1=0; i1<sys->npart; i1++){
    
    if ( ptr[i1].type != types[0] && ptr[i1].type != types[1] )
      continue;

    int moli_i1 = ptr[i1].molindex;

    n = 0;
    while (1) {

      i2 = ptr[i1].neighb[n];
      if ( i2 == -1 ) break; 

      if ( (ptr[i1].type == types[0] && ptr[i2].type == types[1]) || 
         (ptr[i1].type == types[1] && ptr[i2].type == types[0]) ){
	r2 = 0.0;
	for ( k=0; k<3; k++ ){
	  r[k] = ptr[i1].x[k] - ptr[i2].x[k];
	  sep_Wrap( r[k], sys->length[k] );
	  r2 += r[k]*r[k];
	}

	if ( r2 < cf2 ){

	  ft = (*fun)(r2, 'f');
	  for ( k=0; k<3; k++ ){
	    f[k] = ft*r[k];
	    ptr[i1].f[k] += f[k];
	    ptr[i2].f[k] -= f[k];
	  }
	  retval->epot += (*fun)(r2, 'u'); 

	  for (k=0; k<3; k++)
	    for ( int kk=0; kk<3; kk++ )
	      retval->pot_P[k][kk] += f[k]*r[kk];
							
	  if ( sys->molptr->flag_Fij == 1 ){
	    int moli_i2 = ptr[i2].molindex;
	    if ( moli_i1 != -1 && moli_i2 != -1 ){
	      for ( int k=0; k<3;k++ ){
		sys->molptr->Fij[moli_i1][moli_i2][k] += f[k];
		sys->molptr->Fij[moli_i2][moli_i1][k] -= f[k];	
								
		sys->molptr->Fiajb[i1][moli_i2][k] += f[k];
		sys->molptr->Fiajb[i2][moli_i1][k] -= f[k];	
	      }
	    }	
	  }
	  
	}
      }
      n++;
    }
  }
}



void sep_force_pairs(seppart *ptr, const char *types, double cf,
		     double (*fun)(double, char), sepsys *sys, 
		     sepret *retval, const unsigned opt){

  if ( cf > sys->cf )
    sep_error("cutoff for an interaction cannot be larger than maximum cutoff");
  
  sepsys tmp;
  
  if ( sys->neighb_update == SEP_BRUTE ) {
    if (  sys->omp_flag == 1 ){
      sep_warning("omp flag set, SEP_BRUTE does not support threads.");
      sep_warning("Resetting omp flag");
      sys->omp_flag = 0;
    }
    sep_force_pair_brute(ptr, types, cf, fun, sys, retval, opt);    
  }
  else {
    if ( sys->neighb_flag == 1 ){
      
      if ( sys->neighb_update == SEP_NEIGHBLIST ){
	sep_make_neighblist(ptr, sys, opt);
      }
      else if ( sys->neighb_update == SEP_LLIST_NEIGHBLIST ){
	// copying needs to be done to call sep_neighb - should be changed
	memcpy(&tmp, sys, sizeof(sepsys));
	
	if ( opt == SEP_ALL )
	  sep_neighb(ptr, tmp);
	else if ( opt == SEP_NEIGHB_EXCL_BONDED )
	  sep_neighb_nonbonded(ptr, tmp);
	else if ( opt == SEP_NEIGHB_EXCL_SAME_MOL )
	  sep_neighb_excl_same_mol(ptr, tmp);
      }

      sys->neighb_flag = 0;
    }
        
    if ( sys->omp_flag == 1 )
      sep_force_pair_neighb_omp(ptr, types, cf, fun, sys, retval);
    else 
      sep_force_pair_neighb(ptr, types, cf, fun, sys, retval);

  }

}

// Making a neighb list
void sep_make_neighblist(seppart *ptr, sepsys *sys, const unsigned opt){
  int n, m, k, index;
  double r[3], r2, cf2; 
  
  cf2 = sep_Sq(sys->cf + sys->skin); 

  for ( n=0; n<sys->npart; n++ )
    for ( m=0; m<SEP_NEIGHB; m++ ) ptr[n].neighb[m] = -1;

  for ( n=0; n<sys->npart-1; n++ ){
    index = 0;
    for ( m=n+1; m<sys->npart; m++ ){
      
      if ( opt == SEP_NEIGHB_EXCL_BONDED && 
           sep_bonded_direct(ptr, n, m) == 1 ){
        continue;
      }
      else if ( opt == SEP_NEIGHB_EXCL_SAME_MOL && 
                ptr[n].molindex == ptr[m].molindex &&
                ptr[n].molindex != -1 ){
        continue;
      }
      
      r2 = 0.0;
      for ( k=0; k<3; k++ ){
        r[k]  = ptr[n].x[k]-ptr[m].x[k];
        sep_Wrap( r[k], sys->length[k] );
        r2   += r[k]*r[k];
      }

      if ( r2 < cf2 ){ 
        ptr[n].neighb[index] = m;
        index++;
      }
      
    }
  }

}

// LLIST + neigbour list
void sep_neighb_excl_same_mol(seppart *ptr, sepsys sys){
  int *list;

  list = sep_allocate_celllist(sys);
  sep_make_celllist(ptr, list, sys);
  sep_make_neighblist_from_llist_excl_same_mol(ptr, SEP_NEIGHB, list, sys);
  free(list);

}


void sep_neighb(seppart *ptr, sepsys sys){
  int *list;

  list = sep_allocate_celllist(sys);
  sep_make_celllist(ptr, list, sys);
  sep_make_neighblist_from_llist(ptr, SEP_NEIGHB, list, sys);

  free(list);

}


void sep_neighb_nonbonded(seppart *ptr, sepsys sys){
  int *list;

  list = sep_allocate_celllist(sys);
  sep_make_celllist(ptr, list, sys);
  sep_make_neighblist_from_llist_nonbonded(ptr, SEP_NEIGHB, list, sys);

  free(list);

}


int *sep_allocate_celllist(sepsys sys){
  int *ptr;

  ptr = sep_vector_int(sys.npart+
		       sys.nsubbox[0]*sys.nsubbox[1]*sys.nsubbox[2]);

  if ( ptr == NULL )
    sep_error("%s at %d: Memory allocation error");

  return ptr;
}


void sep_make_celllist(seppart *ptr, int *list, sepsys sys){
  int n, i, nsubbox3, nsubbox2, length;
  
  nsubbox2 = sys.nsubbox[0]*sys.nsubbox[1];
  nsubbox3 = nsubbox2*sys.nsubbox[2];
  length = sys.npart + nsubbox3;

  for ( n=0; n<length; n++ ) list[n] = -1;

  for (n=0; n<sys.npart; n++){
    i = (int)(ptr[n].x[0]/sys.lsubbox[0]) + 
      (int)(ptr[n].x[1]/sys.lsubbox[1])*sys.nsubbox[0] +
      (int)(ptr[n].x[2]/sys.lsubbox[2])*nsubbox2;

    if ( i > length - 1 ){
      sep_error("%s at %d: Index larger than array length\n", 
		__func__, __LINE__);
    }

    list[n+nsubbox3] = list[i];
    list[i] = n;  
  }

}


void sep_make_neighblist_from_llist(seppart *ptr,  int nneighb, 
				    int *list, sepsys sys){
  double r[3], r2, cf2;
  int j1, j2, m1, m1X, m1Y, m1Z, m2, m2X, m2Y, m2Z,
    n, k, offset, nsubbox3, nsubbox2;
  static int iofX[] = {0,1,1,0,-1,0,1,1,0,-1,-1,-1,0,1}; 
  static int  iofY[] = {0,0,1,1,1,0,0,1,1,1,0,-1,-1,-1};
  static int  iofZ[] = {0,0,0,0,0,1,1,1,1,1,1,1,1,1}; 
  int *index;

  index = sep_vector_int(sys.npart);  
  nsubbox2 = sys.nsubbox[1]*sys.nsubbox[0]; 
  nsubbox3 = sys.nsubbox[2]*nsubbox2;
  cf2 = sep_Sq(sys.cf + sys.skin);
    
  for (n=0; n<sys.npart; n++)
    for (k=0; k<nneighb; k++)  ptr[n].neighb[k] = -1;
  
    
  for (m1Z = 0; m1Z < sys.nsubbox[2]; m1Z++){
    for (m1Y = 0; m1Y < sys.nsubbox[1]; m1Y++) {
      for (m1X = 0; m1X < sys.nsubbox[0]; m1X++) {
	m1 = m1Z*nsubbox2 + m1Y*sys.nsubbox[0] + m1X;
	for (offset = 0; offset < 14; offset++) {

	  m2X = m1X + iofX[offset];    
	  if (m2X == sys.nsubbox[0] ) m2X = 0;    
	  else if (m2X == -1 ) m2X = sys.nsubbox[0]-1;    
    

	  m2Y = m1Y + iofY[offset];    
	  if (m2Y == sys.nsubbox[1] )  m2Y = 0;    
	  else if (m2Y == -1 ) m2Y = sys.nsubbox[1]-1;    
	  
	  m2Z = m1Z + iofZ[offset];    
	  if ( m2Z == sys.nsubbox[2] ) m2Z = 0;    

	  j1 = list[m1];
	  m2 = m2Z*nsubbox2 + m2Y*sys.nsubbox[0] + m2X;
	  while (j1 != -1) {
	    j2 = list[m2];
	    while (j2 !=-1) {
	      if (m1 != m2 || j2 < j1) {
		r2 = 0.0;
		for (k = 0; k<3; k++){
		  r[k] = ptr[j1].x[k] - ptr[j2].x[k];
		  sep_Wrap( r[k], sys.length[k] );
		  r2 += r[k]*r[k];
		}

		if (r2 < cf2) {
		  ptr[j1].neighb[index[j1]] = j2;
		  index[j1]++;

		  if ( index[j1] > nneighb - 1 ){
		    sep_error("%s at %d: Too many neighbours\n",
			      __func__, __LINE__);
		  }
		  if ( sys.omp_flag == 1 ){
		    ptr[j2].neighb[index[j2]] = j1;
		    index[j2]++;
		  }

		}
	      }
	      j2 = list[j2+nsubbox3];
	    }
	    j1 = list[j1+nsubbox3];  
	  } } } } }

  free(index);
}


void sep_make_neighblist_from_llist_nonbonded(seppart *ptr, int nneighb, 
					      int *list, sepsys sys){
  double r[3], r2, cf2;
  int j1, j2, m1, m1X, m1Y, m1Z, m2, m2X, m2Y, m2Z,
    n, k, offset, nsubbox3, nsubbox2;
  static int iofX[] = {0,1,1,0,-1,0,1,1,0,-1,-1,-1,0,1}; 
  static int iofY[] = {0,0,1,1,1,0,0,1,1,1,0,-1,-1,-1};
  static int iofZ[] = {0,0,0,0,0,1,1,1,1,1,1,1,1,1}; 
  int *index, **bond;
  
  bond = sep_matrix_int(sys.npart, SEP_BOND);
  sep_get_bond_alltoall(ptr, bond, sys.npart);
  
  index = sep_vector_int(sys.npart);  
  nsubbox2 = sys.nsubbox[1]*sys.nsubbox[0]; 
  nsubbox3 = sys.nsubbox[2]*nsubbox2;
  cf2 = sep_Sq(sys.cf + sys.skin);
    
  for (n=0; n<sys.npart; n++){
    for (k=0; k<nneighb; k++)
      ptr[n].neighb[k] = -1;
  }
    
  for (m1Z = 0; m1Z < sys.nsubbox[2]; m1Z++){
    for (m1Y = 0; m1Y < sys.nsubbox[1]; m1Y++) {
      for (m1X = 0; m1X < sys.nsubbox[0]; m1X++) {
	m1 = m1Z*nsubbox2 + m1Y*sys.nsubbox[0] + m1X;
	for (offset = 0; offset < 14; offset++) {
	 
	  m2X = m1X + iofX[offset];    
	  if (m2X == sys.nsubbox[0] ) m2X = 0;    
	  else if (m2X == -1 ) m2X = sys.nsubbox[0]-1;    
    

	  m2Y = m1Y + iofY[offset];    
	  if (m2Y == sys.nsubbox[1] )  m2Y = 0;    
	  else if (m2Y == -1 ) m2Y = sys.nsubbox[1]-1;    
	  
	  m2Z = m1Z + iofZ[offset];    
	  if ( m2Z == sys.nsubbox[2] ) m2Z = 0;    
	  

	  j1 = list[m1];
	  m2 = m2Z*nsubbox2 + m2Y*sys.nsubbox[0] + m2X;
	  while ( j1 != -1 ) {
	    j2 = list[m2];
            while ( j2 != -1 ) {
	      if ( m1 != m2 || j2 < j1 ) {
		r2 = 0.0;
		for ( k = 0; k<3; k++ ){
		  r[k] = ptr[j1].x[k] - ptr[j2].x[k];
		  sep_Wrap( r[k], sys.length[k] );
		  r2 += r[k]*r[k];
		}
		if ( r2 < cf2 && sep_bonded(bond, j1, j2) == 0 ) {
		  ptr[j1].neighb[index[j1]] = j2;
		  index[j1]++;
		  if ( index[j1] > nneighb - 1 ){
		    sep_error("%s at %d: Too many neighbours\n",
			      __func__, __LINE__);
		  }
		  if ( sys.omp_flag == 1 ){
		    ptr[j2].neighb[index[j2]] = j1;
		    index[j2]++;
		  }
		}	
	      }
	      j2 = list[j2+nsubbox3];
	    }
	    j1 = list[j1+nsubbox3];  
	  }
        } } } }

  free(index);
  sep_free_matrix_int(bond, sys.npart);
  
}


void sep_make_neighblist_from_llist_excl_same_mol(seppart *ptr, int nneighb, 
				       int *list, sepsys sys){
  double r[3], r2, cf2;
  int j1, j2, m1, m1X, m1Y, m1Z, m2, m2X, m2Y, m2Z,
    n, k, offset, nsubbox3, nsubbox2;
  static int iofX[] = {0,1,1,0,-1,0,1,1,0,-1,-1,-1,0,1}; 
  static int iofY[] = {0,0,1,1,1,0,0,1,1,1,0,-1,-1,-1};
  static int iofZ[] = {0,0,0,0,0,1,1,1,1,1,1,1,1,1}; 
  int *index;
  
  index = sep_vector_int(sys.npart);  
  nsubbox2 = sys.nsubbox[1]*sys.nsubbox[0]; 
  nsubbox3 = sys.nsubbox[2]*nsubbox2;
  cf2 = sep_Sq(sys.cf + sys.skin);
    
  for (n=0; n<sys.npart; n++){
    for (k=0; k<nneighb; k++)
      ptr[n].neighb[k] = -1;
  }
    
  for (m1Z = 0; m1Z < sys.nsubbox[2]; m1Z++){
    for (m1Y = 0; m1Y < sys.nsubbox[1]; m1Y++) {
      for (m1X = 0; m1X < sys.nsubbox[0]; m1X++) {
	m1 = m1Z*nsubbox2 + m1Y*sys.nsubbox[0] + m1X;
	for (offset = 0; offset < 14; offset++) {
	 
	  m2X = m1X + iofX[offset];    
	  if (m2X == sys.nsubbox[0] ) m2X = 0;    
	  else if (m2X == -1 ) m2X = sys.nsubbox[0]-1;    
    
	  m2Y = m1Y + iofY[offset];    
	  if (m2Y == sys.nsubbox[1] )  m2Y = 0;    
	  else if (m2Y == -1 ) m2Y = sys.nsubbox[1]-1;    
	  
	  m2Z = m1Z + iofZ[offset];    
	  if ( m2Z == sys.nsubbox[2] ) m2Z = 0;    

	  j1 = list[m1];
	  m2 = m2Z*nsubbox2 + m2Y*sys.nsubbox[0] + m2X;
	  while ( j1 != -1 ) {
	    j2 = list[m2];
            while ( j2 != -1 ) {
	      if ( m1 != m2 || j2 < j1 ) {
		
		if ( ptr[j1].molindex == -1 ||  // Not in a molecule 
		     ptr[j1].molindex != ptr[j2].molindex ) {

		  r2 = 0.0;
		  for ( k = 0; k<3; k++ ){
		    r[k] = ptr[j1].x[k] - ptr[j2].x[k];
		    sep_Wrap( r[k], sys.length[k] );
		    r2 += r[k]*r[k];
		  }
		  if ( r2 < cf2  ) {
		    ptr[j1].neighb[index[j1]] = j2;
		    index[j1]++;
		    if ( index[j1] > nneighb - 1 ){
		      sep_error("%s at %d: Too many neighbours\n",
				__func__, __LINE__);
		    }
		    if ( sys.omp_flag == 1 ){
		      ptr[j2].neighb[index[j2]] = j1;
		      index[j2]++;
		    }
		  }

		}

	      }
	      j2 = list[j2+nsubbox3];
	    }
	    j1 = list[j1+nsubbox3];  
	  }
        } } } }

  free(index);
   
}


void sep_get_bond_alltoall(seppart *ptr, int **bond, int npart){
  int i1, i2, k, *index;
  
  index = sep_vector_int(npart);
  sep_matrix_int_set(bond, npart, SEP_BOND, -1);
  for ( i1=0; i1<npart; i1++ ){
    for ( k=0; k<SEP_BOND; k++ ){
      i2 = ptr[i1].bond[k];
      if ( i2 != -1 ){
        bond[i1][index[i1]] = i2;
        bond[i2][index[i2]] = i1;
        index[i1]++;
        index[i2]++;
        if ( index[i1] > SEP_BOND || index[i2] > SEP_BOND )
          sep_error("%s at %d: Index exceeds the number of bonds",
		    __func__, __LINE__);
      }
    }
  }
  
  free(index);        
}

unsigned int sep_bonded(int **bond, int j1, int j2) {
  int k;
  
  for ( k=0; k<SEP_BOND; k++ ) {   
    if ( bond[j1][k] == j2 || bond[j2][k] == j1 ) 
      return 1;
  }
  
  return 0;
}

unsigned int sep_bonded_direct(seppart *ptr, int j1, int j2) {

  for ( int k=0; k<SEP_BOND; k++ ) {   
    if ( ptr[j1].bond[k] == j2 || ptr[j2].bond[k] == j1 ) 
      return 1;
  }
  
  return 0;
}


void sep_force_pair_neighb_omp(seppart *ptr, const char *types, double cf,
			       double (*fun)(double, char), sepsys *sys, 
			       sepret *retval) {
  int i1, i2, n, k;
	double r[3], r2, ft, f[3], epot;
  const double cf2 = cf*cf;

  epot = 0.0;

#pragma omp parallel for			\
  private(n, k, i2, r, r2, ft, f)		\
  reduction(+:epot) 
  for ( i1=0; i1<sys->npart; i1++ ){

    if ( ptr[i1].type != types[0] && ptr[i1].type != types[1] )
      continue;

    n = 0;
    while (1){
      i2 = ptr[i1].neighb[n];
      if ( i2 == -1 ) break; 

      if ( (ptr[i1].type == types[0] && ptr[i2].type == types[1]) || 
	   (ptr[i1].type == types[1] && ptr[i2].type == types[0]) ){
				
	for ( k=0; k<3; k++ ){
	  r[k] = ptr[i1].x[k] - ptr[i2].x[k];
	  sep_Wrap( r[k], sys->length[k] );
	}
	
	r2 = r[0]*r[0] + r[1]*r[1] + r[2]*r[2];
	
	if ( r2 < cf2 ){
	  ft = (*fun)(r2, 'f');
	  for ( k=0; k<3; k++ ){
	    f[k] = ft*r[k];
	    ptr[i1].f[k] += f[k];
	  }
	  epot += (*fun)(r2, 'u');
	}
	
      }
      n++;
    }
  }
  
  retval->epot += 0.5*epot;
}


//////////////////////////////////////////////////////////////////////


void sep_add_interaction(sepsys *sys, char types[], double cf, 
			 double (*fun)(double,char)){
  
  if( cf>sys->cf ) 
    sep_error("%s at %d cutoff of an interaction cannot be larger than maximum cutoff", 
	      __func__, __LINE__);
  
  sys->numinteractions += 1;
  int i =  sys->numinteractions;

  if ( i > 1 ){ // did malloc in sep_setup_sys
    sys->interactionptr = 
      realloc(sys->interactionptr, i*sizeof(sepinteraction));
 
    if ( sys->interactionptr == NULL )
      sep_error("%s at %d: Failure allocation memory", __func__, __LINE__);
  }

  sys->interactionptr[i-1].types[0]=types[0];
  sys->interactionptr[i-1].types[1]=types[1];

  sys->interactionptr[i-1].cf=cf;
  sys->interactionptr[i-1].fun=fun;

}

void _sep_force_pairs(seppart *ptr, sepsys *sys,  sepret *retval, unsigned n, 
		      const unsigned opt){

  sepsys tmp;

  if ( sys->neighb_update == 0 ) {
    _sep_force_pair_brute(ptr, sys, retval, opt);    
  }
  else {
 
    if ( sys->neighb_flag != 0 && (n%sys->neighb_flag==0 || sys->neighb_flag==1) ){
    
      // copying needs to be done to call sep_neighb - should be changed
      memcpy(&tmp, sys, sizeof(sepsys));

      if ( opt == SEP_NEIGHB_ALL )
	sep_neighb(ptr, tmp);
      else if ( opt == SEP_NEIGHB_EXCL_BONDED )
	sep_neighb_nonbonded(ptr, tmp);
      else if ( opt == SEP_NEIGHB_EXCL_SAME_MOL )
	sep_neighb_excl_same_mol(ptr, tmp);

      sys->neighb_flag = 0;
    }
    else if ( (n+1)%sys->neighb_update==0 ){
      sys->neighb_flag = n+1;
    }

    if ( sys->omp_flag == 1 )
      _sep_force_pair_neighb_omp(ptr,  sys, retval);
    else  
    _sep_force_pair_neighb(ptr, sys, retval);
  
  }

}


void _sep_force_pair_brute(seppart *ptr, sepsys *sys, 
			  sepret *retval, const int opt){
  int n,m,k;
  double r[3], r2, ft, f[3], cf2; 
  int i3;
  double cf;
  double (*fun) (double,char);
  char types[2];
  

  for (n=0; n<sys->npart-1; n++){
		
    int i = ptr[n].molindex;
					
    for (m=n+1; m<sys->npart; m++){

      if ( opt == SEP_NEIGHB_EXCL_BONDED && sep_bonded_direct(ptr, n, m) == 1 ){
	continue;
      }
      else if ( opt == SEP_NEIGHB_EXCL_SAME_MOL && 
		ptr[n].molindex == ptr[m].molindex &&
		ptr[n].molindex != -1 )
	continue;
      
      for ( i3=0; i3<sys->numinteractions; i3++ ){
	types[0] = sys->interactionptr[i3].types[0];
	types[1] = sys->interactionptr[i3].types[1];
	
	cf  = sys->interactionptr[i3].cf;
	fun = sys->interactionptr[i3].fun;
	cf2 = cf*cf;

	if ( (ptr[n].type == types[0] && ptr[m].type == types[1]) || 
	     (ptr[n].type == types[1] && ptr[m].type == types[0]) ){
	  
	  r2 = 0.0;
	  for (k=0; k<3; k++){
	    r[k]  = ptr[n].x[k]-ptr[m].x[k];
	    sep_Wrap( r[k], sys->length[k] );
	    r2   += r[k]*r[k];
	  }

	  if (r2 < cf2){ 

	    // Force between particles	
	    ft = (*fun)(r2, 'f'); 
	    
	    for (k=0; k<3; k++){
	      f[k] = ft*r[k];
	      
	      ptr[n].f[k] += f[k];
	      ptr[m].f[k] -= f[k];
	    }
	    
	    // Energy 
	    retval->epot += (*fun)(r2, 'u');
	    
	    // Config. part of the stress/pressure tensor 
	    for (k=0; k<3; k++)
	      for ( int kk=0; kk<3; kk++ ) retval->pot_P[k][kk] += f[k]*r[kk];
	    
	    // Force between molecules
	    if ( sys->molptr->flag_Fij == 1 ){
	      int j = ptr[m].molindex;
	      if ( i != -1 && j != -1 ){
		
		for ( k=0; k<3;k++ ){
		  // pressure tensor	
		  sys->molptr->Fij[i][j][k] += f[k];
		  sys->molptr->Fij[j][i][k] -= f[k];
		  // Couple tensor
		  sys->molptr->Fiajb[n][j][k] += f[k];
		  sys->molptr->Fiajb[m][i][k] -= f[k];	
		}
	      }	
	    }
	    
	  } // End of if ( r2 < rcut )
	}
      }
    }
  }
  
}

void _sep_force_pair_neighb(seppart *ptr,  sepsys *sys, 
			    sepret *retval) {
  int i1, i2, i3, n, k;
  double r2, ft, f[3], r[3];
  double cf2, cf;
  double (*fun) (double,char);
  char types[2];
  
  for (i1=0; i1<sys->npart; i1++){
   for (i3=0; i3<sys->numinteractions; i3++){
    types[0]=sys->interactionptr[i3].types[0];
    types[1]=sys->interactionptr[i3].types[1];
    
    cf  = sys->interactionptr[i3].cf;
    fun = sys->interactionptr[i3].fun;
    
    if ( ptr[i1].type != types[0] && ptr[i1].type != types[1] )
      continue;

    int moli_i1 = ptr[i1].molindex;

    n = 0;
    while (1){
      i2 = ptr[i1].neighb[n];
      if ( i2 == -1 ) break; 

      if ( (ptr[i1].type == types[0] && ptr[i2].type == types[1]) || 
	   (ptr[i1].type == types[1] && ptr[i2].type == types[0]) ){
	r2 = 0.0;
	for ( k=0; k<3; k++ ){
	  r[k] = ptr[i1].x[k] - ptr[i2].x[k];
	  sep_Wrap( r[k], sys->length[k] );
	  r2 += r[k]*r[k];
	}
        cf2=cf*cf;
	if ( r2 < cf2 ){
	  ft = (*fun)(r2, 'f');
	  for ( k=0; k<3; k++ ){
	    f[k] = ft*r[k];
	    ptr[i1].f[k] += f[k];
	    ptr[i2].f[k] -= f[k];
	  }
	  retval->epot += (*fun)(r2, 'u'); 

	  for (k=0; k<3; k++)
	    for ( int kk=0; kk<3; kk++ )
	      retval->pot_P[k][kk] += f[k]*r[kk];
							
	  if ( sys->molptr->flag_Fij == 1 ){
	    int moli_i2 = ptr[i2].molindex;
	    if ( moli_i1 != -1 && moli_i2 != -1 ){
	      for ( int k=0; k<3;k++ ){
		sys->molptr->Fij[moli_i1][moli_i2][k] += f[k];
		sys->molptr->Fij[moli_i2][moli_i1][k] -= f[k];	
								
		sys->molptr->Fiajb[i1][moli_i2][k] += f[k];
		sys->molptr->Fiajb[i2][moli_i1][k] -= f[k];	

	      }
	    }	
	  }

	}
      }
      n++;
    }
   }
  }

}


void _sep_force_pair_neighb_omp(seppart *ptr, sepsys *sys,  sepret *retval) {
  int i1, i2, i3, n, k;
  double r[3], r2, ft, f[3], epot, cf2;
  char t0, t1;

  epot = 0.0;
#pragma omp parallel for				\
  private(n, k, i2, i3, cf2, r, r2, ft, f, t0, t1)	\
  reduction(+:epot) 
  for ( i1=0; i1<sys->npart; i1++ ){

    n = 0;
    while (1){
      i2 = ptr[i1].neighb[n];
      if ( i2 == -1 ) break; 

      for (i3=0; i3<sys->numinteractions; i3++){
	t0=sys->interactionptr[i3].types[0];
	t1=sys->interactionptr[i3].types[1];
    
	cf2  = sep_Sq(sys->interactionptr[i3].cf);
	
	if ( (ptr[i1].type == t0 && ptr[i2].type == t1) || 
	     (ptr[i1].type == t1 && ptr[i2].type == t0) ){
				
	  for ( k=0; k<3; k++ ){
	    r[k] = ptr[i1].x[k] - ptr[i2].x[k];
	    sep_Wrap( r[k], sys->length[k] );
	  }
	  
	  r2 = r[0]*r[0] + r[1]*r[1] + r[2]*r[2];
	  
	  if ( r2 < cf2 ){
	    ft = sys->interactionptr[i3].fun(r2, 'f');
	    for ( k=0; k<3; k++ ){
	      f[k] = ft*r[k];
	      ptr[i1].f[k] += f[k];
	    }
	    epot += sys->interactionptr[i3].fun(r2, 'u');
	  }
	  
	}
      }
      n++;
    }
  }
  
  retval->epot += 0.5*epot;
}

	
